#define THEOR_BI_GLOBAL_VARIABLES
#include "theor_bispectrum_5.h"

#define TILT_ADI 0.963
#define TILT_ISO 0.963
#define TILT_ISO_STRING "0.963_indep_"
/*
0.963
0.5, 1.5
0.7, 1.2
0.0, 2.0
0.9, 1.0
0.3, 1.7
3.0, 4.0, etc.
*/
#define AMPL 2.43e-9    // for WMAP7-only
// define AMPL 2.441e-9    // for WMAP7 + H0 + BAO
#define PIVOT 0.002
#define REL_AMPL_AI 0.0016
#define REL_AMPL_ISO 0.014
#define NORM 2.*M_PI*M_PI * 9./25. * AMPL
#define REALRUN 3 /* 0 = print some output for testing, no files are written
		     1 = do power spectrum only
		     2 = power spectrum + bispectrum
		     3 = independant case */
#define WRITE_KINTVARS 0 // For testing purposes (1 = yes, 0 = no)

#define NUM_SHAPE 3 // Number of bispectrum shapes implemented
int do_shape = 1;   /* Which bispectrum shapes to run
		       0 = all
		       1 = local only
		       2 = equilateral only
		       4 = orthogonal only
		       Take sums for combinations
		       Note: unresolved extra-galactic point sources have
		       a constant bispectrum which is taken into account
		       directly in bispec_calc, no need to include it here. */

#define NUM_ISO 5 // Number of isocurvature components implemented
int do_iso[NUM_SHAPE] = {3, 1, 1};
// For the moment isocurvature is not supported for non-local shapes
/* Which isocurvature components to run for each shape
   0 = all
   1 = adiabatic only
   2 = CDM isocurvature only
   4 = neutrino density isocurvature only
   8 = neutrino velocity isocurvature only
   16 = baryon isocurvature only
   Take sums for combinations */

#define NUM_POL 2 // Number of polarization components implemented
int do_pol[NUM_SHAPE] = {0, 0, 0};
/* Which polarization components to run for each shape
   0 = all
   1 = temperature only
   2 = E-polarization only
   Take sums for combinations */

const double transfer_corr[NUM_ISO] = {5./3., 5., 5., 5., 5.};
// The correction factor to convert CAMB's transfer function to Phi as
// used in this code.
const double rel_ampl_list[NUM_ISO*(NUM_ISO+1)/2] =
//  {1., 0., 0., 0., 0., 1./9., 0., 0., 0., 1./9., 0., 0., 1./9., 0., 1./9.};
  {1., 1./3., 1./3., 1./3., 1./3., 1./9., 1./9., 1./9., 1./9.,
   1./9., 1./9., 1./9., 1./9., 1./9., 1./9.};
 //{1.,  REL_AMPL_AI,  REL_AMPL_AI,  REL_AMPL_AI,  REL_AMPL_AI,
 //      REL_AMPL_ISO, REL_AMPL_ISO, REL_AMPL_ISO, REL_AMPL_ISO,
 //                    REL_AMPL_ISO, REL_AMPL_ISO, REL_AMPL_ISO,
 //                                  REL_AMPL_ISO, REL_AMPL_ISO,
 //                                                REL_AMPL_ISO};

/* Relative amplitude of isocurvature and cross correlated
   power spectra with respect to adiabatic one.
   Order: adi-adi, adi-iso, iso-iso;
   or more general: a-a, a-i1, a-i2, a-i3, a-i4, i1-i1, i1-i2, i1-i3, i1-i4,
   i2-i2, i2-i3, i2-i4, i3-i3, i3-i4, i4-i4.
   In general, apart from the first entry which is always 1, these should be
   expressions like AMPL_iso*pow(PIVOT_iso, 1-TILT_iso) divided by
   (AMPL_adi*pow(PIVOT_adi, 1-TILT_adi)).
   In addition there is a factor 1/9 for the isocurvature components, due
   to conversions between Phi and zeta/S. For the adiabatic component
   the correction factor is (3/5)^2, included in NORM. For the isocurvature
   components this factor should be (1/5)^2, i.e. 1/9 times smaller.

   For the moment the cross terms only change the power spectra.
   For the bispectra it is assumed that only the adiabatic primordial
   power spectrum plays a role. However, due to conversions between Phi
   and zeta/S, the code still uses all non-cross primordial power spectra,
   but the relative amplitude of the non-adiabatic ones should be kept
   at 1./9. and the tilt equal to the adiabatic one.
   Once the generalization has been worked out theoretically, it will
   be implemented. It would also probably be a good idea to get rid of these
   Phi/zeta conversion factors at that time, and only work with zeta.
*/
const double red_tilt_list[NUM_ISO*(NUM_ISO+1)/2] =
  {TILT_ADI-1, (TILT_ADI+TILT_ISO)/2.-1, (TILT_ADI+TILT_ISO)/2.-1, (TILT_ADI+TILT_ISO)/2.-1, (TILT_ADI+TILT_ISO)/2.-1,
                             TILT_ISO-1,               TILT_ISO-1,               TILT_ISO-1,               TILT_ISO-1,
                                                       TILT_ISO-1,               TILT_ISO-1,               TILT_ISO-1,
                                                                                 TILT_ISO-1,               TILT_ISO-1,
                                                                                                           TILT_ISO-1};
// Reduced spectral index \tilde{n} = n - 1. Same order as rel_ampl_list.

#define FILENAME_LEN 200
// The filenames below (after construction, where applicable) should be no
// more than filename_len characters.
#define absolute_path "/sps/planck/Users/linc/"
#define cosmology_string "newCamb_3100_"
#define tf_cosmology_string "newCamb_wmap7only_3100_"
#define filename_end ".dat"
#define filename_tf_start "TransferFunction_"
#define filename_bispec_start "Bispectrum_"
#define filename_cl_start "Powerspectrum_"
#define filename_kintvars_start "kintvars_"
#define filename_index_start "Index_"

#define TRIPLET_TYPE 2
/* 0 means that only triplets satisfying the triangle inequality are computed,
   1 means that an additional number of triplets necessary for linear
   interpolation is computed (so that all points inside the triangle can be
   reached with linear interpolation),
   2 is analogous to 1, but for a cubic interpolation scheme,
   3 means that all triplets are computed, not taking into account the
   triangle inequality at all. */
#define PARITY_TYPE 1
/* 0 means that only triplets satisfying the parity condition are computed,
   1 means the parity condition is not taken into account (all triplets
   are included) */
int stepstep = 10; // only matters if stepmax>1 (see function adaptive_step)
int stepmax = 1;   // stepmax=1 means that all points are used
int limin = 0;
int limax = -1;    // Put explicit value, or -1 for last l_index value

#define STRINGLEN 10
const char tria_list[4][STRINGLEN] = {"noextra", "lin", "cub", "notria"};
const char par_list[2][STRINGLEN] = {"par", "nopar"};
const char shape_list[NUM_SHAPE][STRINGLEN] = {"local", "equil", "ortho"};
const char iso_list[NUM_ISO][STRINGLEN] = {"adi", "cdm", "nden", "nvel", "bar"};
const char pol_list[NUM_POL][STRINGLEN] = {"T", "E"};

transferfunc_vars tfvars;
k_integral_vars kintvars;


int triplet_relevant_noextra(int li1, int li2, int li3){
  int l1, l2, l3;

  l1 = tfvars.l_arr[li1];
  l2 = tfvars.l_arr[li2];
  l3 = tfvars.l_arr[li3];
  if (l3 <= l1+l2)
    return 1;
  else
    return 0;
}

int triplet_relevant_lin(int li1, int li2, int li3, int stepstep, int stepmax){
  int l1max, l2max, l3min, li1step, li2step, li3step;

  li1step = adaptive_step(li1,stepstep,stepmax);
  li2step = adaptive_step(li2,stepstep,stepmax);
  li3step = adaptive_step(li3,stepstep,stepmax);
  if (li1+li1step > tfvars.num_l-1)
    l1max = tfvars.l_arr[tfvars.num_l-1];
  else
    l1max = tfvars.l_arr[li1+li1step] - 1;
  if (li2+li2step > tfvars.num_l-1)
    l2max = tfvars.l_arr[tfvars.num_l-1];
  else
    l2max = tfvars.l_arr[li2+li2step] - 1;
  if (li3 < li2+li3step)
    l3min = tfvars.l_arr[li3];
  else
    l3min = tfvars.l_arr[li3-li3step] + 1;
  if (l3min <= l1max+l2max)
    return 1;
  else
    return 0;
}

int triplet_relevant_cub(int li1, int li2, int li3, int stepstep, int stepmax){
  int l1max,l2max,l3min, li1step,li2step,li3step, li1step2,li2step2,li3step2;

  li1step = adaptive_step(li1,stepstep,stepmax);
  li2step = adaptive_step(li2,stepstep,stepmax);
  li3step = adaptive_step(li3,stepstep,stepmax);
  li1step2 = adaptive_step(li1+li1step,stepstep,stepmax);
  li2step2 = adaptive_step(li2+li2step,stepstep,stepmax);
  li3step2 = adaptive_step(li3-li3step,stepstep,stepmax);
  if (li1+li1step+li1step2 > tfvars.num_l-1)
    l1max = tfvars.l_arr[tfvars.num_l-1];
  else
    l1max = tfvars.l_arr[li1+li1step+li1step2] - 1;
  if (li2+li2step+li2step2 > tfvars.num_l-1)
    l2max = tfvars.l_arr[tfvars.num_l-1];
  else
    l2max = tfvars.l_arr[li2+li2step+li2step2] - 1;
  if (li3 < li2+li3step+li3step2)
    l3min = tfvars.l_arr[li3];
  else
    l3min = tfvars.l_arr[li3-li3step-li3step2] + 1;
  if (l3min <= l1max+l2max)
    return 1;
  else
    return 0;
}

int triplet_relevant_wrap(int type, int li1, int li2, int li3, int stepstep, int stepmax){

  if ((PARITY_TYPE == 0) &&
      (((tfvars.l_arr[li1]+tfvars.l_arr[li2]+tfvars.l_arr[li3]) % 2) != 0))
    return 0;
  switch (type) {
  case 0: return triplet_relevant_noextra(li1, li2, li3);
  case 1: return triplet_relevant_lin(li1, li2, li3, stepstep, stepmax);
  case 2: return triplet_relevant_cub(li1, li2, li3, stepstep, stepmax);
  case 3: return 1;
  default:
    fprintf(stderr, "Unsupported triplet selection type %d.\n", type);
    exit(1);
  }
}

void output_kintvars(char iso_str[], char pol_str[]){
// For testing purposes
  FILE *fp;
  int i, j, iso, pol;
  int ni = tfvars.num_iso;
  int np = tfvars.num_pol;
  int nl = kintvars.num_li;
  int nx = kintvars.num_x;
  char filename_out[FILENAME_LEN];

  strcpy(filename_out, absolute_path);
  strcat(filename_out, filename_kintvars_start);
  strcat(filename_out, cosmology_string);
  strcat(filename_out, TILT_ISO_STRING);
  strcat(filename_out, pol_str);
  strcat(filename_out, "_");
  strcat(filename_out, iso_str);
  strcat(filename_out, filename_end);
  fp = fopen(filename_out, "w");
  fprintf(fp, "# num_l = %d\n", nl);
  fprintf(fp, "# num_x = %d\n", nx);
  fprintf(fp, "# num_pol = %d\n", tfvars.num_pol_do);
  fprintf(fp, "# num_iso = %d\n", tfvars.num_iso_do);
  fprintf(fp, "# iso \t pol \t l \t x \t\t alpha \t\t beta \t\t gamma \t\t delta\n");
  for (iso = 0; iso < ni; iso++)
    for (pol = 0; pol < np; pol++)
      if (tfvars.iso_summ[iso] && tfvars.pol_summ[pol])
	for (i = 0; i < nl; i++)
	  for (j = 0; j < nx; j++)
	    if (do_shape > 1 && tfvars.iso_summ_nl[iso]
		&& tfvars.pol_summ_nl[pol])
	      fprintf(fp, "%s \t %s \t %d \t %f \t %e \t %e \t %e \t %e\n",
		      iso_list[iso], pol_list[pol],
		      tfvars.l_arr[kintvars.li_arr[i]], kintvars.x_arr[j],
		      2./M_PI*kintvars.alpha_arr[iso*np*nl*nx+pol*nl*nx+i*nx+j],
		      2./M_PI*NORM*kintvars.beta_adi_arr[iso*np*nl*nx+pol*nl*nx+i*nx+j],
		      2./M_PI*pow(NORM,1./3.)*kintvars.gamma_arr[iso*np*nl*nx+pol*nl*nx+i*nx+j],
		      2./M_PI*pow(NORM,2./3.)*kintvars.delta_arr[iso*np*nl*nx+pol*nl*nx+i*nx+j]);
	    else
	      fprintf(fp, "%s \t %s \t %d \t %f \t %e \t %e \t %e \t %e\n",
		      iso_list[iso], pol_list[pol],
		      tfvars.l_arr[kintvars.li_arr[i]], kintvars.x_arr[j],
		      2./M_PI*kintvars.alpha_arr[iso*np*nl*nx+pol*nl*nx+i*nx+j],
		      2./M_PI*NORM*kintvars.beta_adi_arr[iso*np*nl*nx+pol*nl*nx+i*nx+j],
		      0., 0.);
  fclose(fp);
}

void print_stuff(void){
  int i, j;

  printf("num_l = %d\tnum_p = %d\n", tfvars.num_l, tfvars.num_p);
  printf("num_pol = %d\tnum_iso = %d\tnum_shape = %d\n",
	 tfvars.num_pol, tfvars.num_iso, tfvars.num_shape);
  printf("num_pol_do = %d\tnum_iso_do = %d\tnum_shape_do = %d\n",
	 tfvars.num_pol_do, tfvars.num_iso_do, tfvars.num_shape_do);
  printf("do_pol_summ = %d\tdo_iso_summ = %d\n",
	 tfvars.do_pol_summ, tfvars.do_iso_summ);
  printf("qmin = %g\tqmax = %g\n",tfvars.p_arr[0],tfvars.p_arr[tfvars.num_p-1]);
  printf("l_arr:\n");
  for (i = 0; i < tfvars.num_l; i++)
    printf("%d\t", tfvars.l_arr[i]);
  printf("\n");
  printf("shape_arr:\n");
  for (i = 0; i < tfvars.num_shape; i++)
    printf("%d\t", tfvars.shape_arr[i]);
  printf("\n");
  printf("pol_arr:\n");
  for (i = 0; i < tfvars.num_shape; i++){
    for (j = 0; j < tfvars.num_pol; j++)
      printf("%d\t", tfvars.pol_arr[i][j]);
    printf("\n");
  }
  printf("iso_arr:\n");
  for (i = 0; i < tfvars.num_shape; i++){
    for (j = 0; j < tfvars.num_iso; j++)
      printf("%d\t", tfvars.iso_arr[i][j]);
    printf("\n");
  }
  printf("pol_summ:\n");
  for (i = 0; i < tfvars.num_pol; i++)
    printf("%d\t", tfvars.pol_summ[i]);
  printf("\n");
  printf("pol_summ_nl:\n");
  for (i = 0; i < tfvars.num_pol; i++)
    printf("%d\t", tfvars.pol_summ_nl[i]);
  printf("\n");
  printf("iso_summ:\n");
  for (i = 0; i < tfvars.num_iso; i++)
    printf("%d\t", tfvars.iso_summ[i]);
  printf("\n");
  printf("iso_summ_nl:\n");
  for (i = 0; i < tfvars.num_iso; i++)
    printf("%d\t", tfvars.iso_summ_nl[i]);
  printf("\n");
}

int *which_pol_iso_shape(int incl, int num){
  int i;
  int *arr;

  arr = (int*)malloc(num * sizeof(int));
  for (i = 0; i < num; i++)
    if (incl & ((int)pow(2,i)))
      arr[i] = 1;
    else
      arr[i] = 0;
  return arr;
}

char **make_pol_iso_string(int num, const char list[][STRINGLEN]){
  char **s;
  int i, j;

  s = (char**) malloc((int)pow(2,num) * sizeof(char*));
  for (i = 0; i < (int)pow(2,num); i++){
    s[i] = (char*) malloc(num*STRINGLEN * sizeof(char));
    strcpy(s[i], "");
    for (j = 0; j < num; j++)
      if (i & (int)pow(2,j)){
	strcat(s[i], list[j]);
	strcat(s[i], "+");
      }
    if (i > 0)
      s[i][strlen(s[i])-1] = 0;
  }
  return s;
}

void initialization(void){
  int i, j, k;

  tfvars.pivot = PIVOT;
  tfvars.red_tilt_adi = TILT_ADI - 1;
  tfvars.red_tilt_iso = TILT_ISO - 1;
  tfvars.rel_ampl_ai = REL_AMPL_AI;
  tfvars.rel_ampl_iso = REL_AMPL_ISO;
  tfvars.num_shape = NUM_SHAPE;
  tfvars.num_iso = NUM_ISO;
  tfvars.num_pol = NUM_POL;
  tfvars.num_shape_do = 0;
  tfvars.num_iso_do = 0;
  tfvars.num_pol_do = 0;
  if (do_shape == 0)
    do_shape = (int)pow(2,NUM_SHAPE) - 1;
  tfvars.shape_arr = which_pol_iso_shape(do_shape, NUM_SHAPE);
  tfvars.iso_arr = (int**)malloc(NUM_SHAPE * sizeof(int*));
  tfvars.pol_arr = (int**)malloc(NUM_SHAPE * sizeof(int*));
  tfvars.iso_summ = (int*)malloc(NUM_ISO * sizeof(int));
  tfvars.iso_summ_nl = (int*)malloc(NUM_ISO * sizeof(int));
  tfvars.pol_summ = (int*)malloc(NUM_POL * sizeof(int));
  tfvars.pol_summ_nl = (int*)malloc(NUM_POL * sizeof(int));
  tfvars.tf_corr = (double*)malloc(NUM_ISO * sizeof(double));
  for (j = 0; j < NUM_ISO; j++){
    tfvars.tf_corr[j] = transfer_corr[j];
    tfvars.iso_summ[j] = 0;
    tfvars.iso_summ_nl[j] = 0;
  }
  for (j = 0; j < NUM_POL; j++){
    tfvars.pol_summ[j] = 0;
    tfvars.pol_summ_nl[j] = 0;
  }
  for (i = 0; i < NUM_SHAPE; i++){
    tfvars.num_shape_do += tfvars.shape_arr[i];
    if (do_iso[i] == 0)
      do_iso[i] = (int)pow(2,NUM_ISO) - 1;
    if (do_pol[i] == 0)
      do_pol[i] = (int)pow(2,NUM_POL) - 1;
    tfvars.iso_arr[i] = which_pol_iso_shape(do_iso[i], NUM_ISO);
    tfvars.pol_arr[i] = which_pol_iso_shape(do_pol[i], NUM_POL);
    for (j = 0; j < NUM_ISO; j++){
      tfvars.iso_summ[j] = (tfvars.iso_summ[j] | tfvars.iso_arr[i][j]);
      if (i > 0)
	tfvars.iso_summ_nl[j] = (tfvars.iso_summ_nl[j] | tfvars.iso_arr[i][j]);
    }
    for (j = 0; j < NUM_POL; j++){
      tfvars.pol_summ[j] = (tfvars.pol_summ[j] | tfvars.pol_arr[i][j]);
      if (i > 0)
	tfvars.pol_summ_nl[j] = (tfvars.pol_summ_nl[j] | tfvars.pol_arr[i][j]);
    }
  }
  for (j = 0; j < NUM_ISO; j++)
    tfvars.num_iso_do += tfvars.iso_summ[j];
  for (j = 0; j < NUM_POL; j++)
    tfvars.num_pol_do += tfvars.pol_summ[j];
  tfvars.do_iso_summ = 0;
  tfvars.do_pol_summ = 0;
  for (i = 0; i < NUM_ISO; i++)
    tfvars.do_iso_summ += tfvars.iso_summ[i] * (int)pow(2,i);
  for (i = 0; i < NUM_POL; i++)
    tfvars.do_pol_summ += tfvars.pol_summ[i] * (int)pow(2,i);
  tfvars.red_tilt = (double**) malloc(NUM_ISO * sizeof(double*));
  tfvars.rel_ampl = (double**) malloc(NUM_ISO * sizeof(double*));
  k = 0;
  for (i = 0; i < NUM_ISO; i++){
    tfvars.red_tilt[i] = (double*) malloc(NUM_ISO * sizeof(double));
    tfvars.rel_ampl[i] = (double*) malloc(NUM_ISO * sizeof(double));
    for (j = 0; j < NUM_ISO; j++)
      if (j >= i){
	tfvars.red_tilt[i][j] = red_tilt_list[k];
	tfvars.rel_ampl[i][j] = rel_ampl_list[k];
	k++;
      }
      else{
	tfvars.red_tilt[i][j] = tfvars.red_tilt[j][i];
	tfvars.rel_ampl[i][j] = tfvars.rel_ampl[j][i];
      }
  }
}

void free_tfvars(void){
  int i;

  free(tfvars.shape_arr);
  for (i = 0; i < tfvars.num_shape; i++){
    free(tfvars.iso_arr[i]);
    free(tfvars.pol_arr[i]);
  }
  free(tfvars.iso_arr);
  free(tfvars.pol_arr);
  free(tfvars.iso_summ);
  free(tfvars.pol_summ);
  free(tfvars.iso_summ_nl);
  free(tfvars.pol_summ_nl);
  free(tfvars.tf_corr);
  for (i = 0; i < tfvars.num_iso; i++){
    free(tfvars.rel_ampl[i]);
    free(tfvars.red_tilt[i]);
  }
  free(tfvars.rel_ampl);
  free(tfvars.red_tilt);
  free(tfvars.l_arr);
  free(tfvars.p_arr);
}


int main(void){
  funcparams params;
  int i, j, k, m, lii1, lii2, lii3, count, shape;
  int np, ni, dimi, dimb, index, li3min, first, do_pol_summ, do_iso_summ;
  int l[3], *index_arr;
  double *bispec, *bispectemp, *cl, *cltemp;
  FILE *fp;
  char **filenames_tf, **iso_string, **pol_string;
  char filename_out[FILENAME_LEN], filename_index[FILENAME_LEN],
    filename_temp[FILENAME_LEN];

  initialization();
  iso_string = make_pol_iso_string(NUM_ISO, iso_list);
  pol_string = make_pol_iso_string(NUM_POL, pol_list);
  filenames_tf = (char**) malloc(NUM_ISO * sizeof(char*));
  for (i = 0; i < NUM_ISO; i++){
    filenames_tf[i] = (char*) malloc(FILENAME_LEN * sizeof(char));
    strcpy(filenames_tf[i], absolute_path);
    strcat(filenames_tf[i], filename_tf_start);
    strcat(filenames_tf[i], tf_cosmology_string);
    strcat(filenames_tf[i], pol_string[(int)pow(2,NUM_POL)-1]);
    strcat(filenames_tf[i], "_");
    strcat(filenames_tf[i], iso_string[(int)pow(2,i)]);
    strcat(filenames_tf[i], filename_end);
  }
  init_theor_bispec(filenames_tf);
  if (limax == -1)
    limax = tfvars.num_l - 1;
  params.qmin = tfvars.p_arr[0];
  params.qmax = tfvars.p_arr[tfvars.num_p-1];
  if (!REALRUN)
    print_stuff();

  init_kintvars_x();
  init_kintvars_w();
  init_kintvars_l(limin, limax, stepstep, stepmax, do_shape);
  if (REALRUN > 1 || WRITE_KINTVARS)
    calculate_k_integrals(&params, do_shape);
  if (WRITE_KINTVARS)
    output_kintvars(iso_string[tfvars.do_iso_summ],
		    pol_string[tfvars.do_pol_summ]);

  if (REALRUN){
    np = tfvars.num_pol_do; ni = tfvars.num_iso_do;
    cl = (double*)malloc((limax-limin+1)*np*(np+1)/2*ni*(ni+1)/2*sizeof(double));
    if (cl == NULL){
      fprintf(stderr, "Memory allocation failed.\n");
      exit(1);
    }
    for (lii1 = limin; lii1 <= limax; lii1++){
      params.li[0] = lii1;
      cltemp = compute_cl(&params);
      for (j = 0; j < ni*(ni+1)/2; j++)
	for (i = 0; i < np*(np+1)/2; i++)
	  cl[lii1-limin + i*(limax-limin+1) + j*(limax-limin+1)*np*(np+1)/2]
	    = cltemp[i + j*np*(np+1)/2];
      free(cltemp);
    }
    do_pol_summ = 0;
    for (i = 0; i < NUM_POL; i++)
      if (tfvars.pol_summ[i])
	do_pol_summ += (int)pow(2,i);
    do_iso_summ = 0;
    for (i = 0; i < NUM_ISO; i++)
      if (tfvars.iso_summ[i])
	do_iso_summ += (int)pow(2,i);
    strcpy(filename_out, absolute_path);
    strcat(filename_out, filename_cl_start);
    strcat(filename_out, cosmology_string);
    strcat(filename_out, TILT_ISO_STRING);
    strcat(filename_out, pol_string[tfvars.do_pol_summ]);
    strcat(filename_out, "_");
    strcat(filename_out, iso_string[tfvars.do_iso_summ]);
    strcat(filename_out, filename_end);
    fp = fopen(filename_out, "w");
    fprintf(fp, "# Header: num_pol %i\tnum_iso %i\tincl_pol %i\tincl_iso %i\n",
	    np, ni, do_pol_summ, do_iso_summ);
    fprintf(fp, "# Power spectra multiplied by l(l+1)\n");
    fprintf(fp, "# l\t");
    for (k = 0; k < NUM_ISO; k++)
      for (m = k; m < NUM_ISO; m++)
	for (i = 0; i < NUM_POL; i++)
	  for (j = i; j < NUM_POL; j++)
	    if (tfvars.iso_summ[k] && tfvars.iso_summ[m]
		&& tfvars.pol_summ[i] && tfvars.pol_summ[j])
	      fprintf(fp, "%s-%s|%s%s\t\t",
		      iso_list[k], iso_list[m], pol_list[i], pol_list[j]);
    fprintf(fp, "\n#\n");
    for (lii1 = limin; lii1 <= limax; lii1++){
      l[0] = tfvars.l_arr[lii1];
      fprintf(fp, "%d\t", l[0]);
      for (i = 0; i < np*(np+1)/2 * ni*(ni+1)/2; i++)
	fprintf(fp, "%.12e\t",
		cl[lii1-limin + i*(limax-limin+1)] * l[0]*(l[0]+1) * NORM);
      fprintf(fp, "\n");
    }
    fclose(fp);
    printf(" *** Successfully create power spectrum. *** \n");
  }

  count = 0;
  for (lii1 = 0; lii1 < kintvars.num_li; lii1++)
    for (lii2 = lii1; lii2 < kintvars.num_li; lii2++)
      for (lii3 = lii2; lii3 < kintvars.num_li; lii3++)
	if (triplet_relevant_wrap(TRIPLET_TYPE, kintvars.li_arr[lii1],
				  kintvars.li_arr[lii2], kintvars.li_arr[lii3],
				  stepstep, stepmax)){
	  count++;
	  //if (!REALRUN){
	  //  l[0] = tfvars.l_arr[kintvars.li_arr[lii1]];
	  //  l[1] = tfvars.l_arr[kintvars.li_arr[lii2]];
	  //  l[2] = tfvars.l_arr[kintvars.li_arr[lii3]];
	  //  printf("l1 = %i\tl2 = %i\tl3 = %i\n", l[0], l[1], l[2]);
	  //}
	}
  if (!REALRUN)
    printf("Total number of triplets: %d\n", count);

  if (REALRUN > 1)
    for (shape = 0; shape < NUM_SHAPE; shape++)
      if (tfvars.shape_arr[shape]){
	strcpy(filename_out, absolute_path);
	strcat(filename_out, filename_bispec_start);
	strcpy(filename_index, absolute_path);
	strcat(filename_index, filename_index_start);
	strcpy(filename_temp, cosmology_string);
	strcat(filename_temp, TILT_ISO_STRING);
	strcat(filename_temp, tria_list[TRIPLET_TYPE]);
	strcat(filename_temp, "_");
	strcat(filename_temp, par_list[PARITY_TYPE]);
	strcat(filename_temp, "_");
	strcat(filename_temp, pol_string[do_pol[shape]]);
	strcat(filename_temp, "_");
	strcat(filename_temp, iso_string[do_iso[shape]]);
	strcat(filename_temp, "_");
	strcat(filename_temp, shape_list[shape]);
	strcat(filename_temp, filename_end);
	strcat(filename_out, filename_temp);
	strcat(filename_index, filename_temp);

	np = 0; ni = 0;
	for (i = 0; i < NUM_POL; i++)
	  np += tfvars.pol_arr[shape][i];
	for (i = 0; i < NUM_ISO; i++)
	  ni += tfvars.iso_arr[shape][i];
	dimi = 3*kintvars.num_li*(kintvars.num_li+1)/2 + kintvars.num_li + 10;
	index_arr = (int*)malloc(dimi * sizeof(int));
	index_arr[0] = tfvars.l_arr[kintvars.li_arr[0]];
	index_arr[1] = tfvars.l_arr[kintvars.li_arr[kintvars.num_li-1]];
	index_arr[2] = kintvars.num_li;
	index_arr[3] = count;
	index_arr[4] = np;
	index_arr[5] = ni;
	index_arr[6] = do_pol[shape];
	index_arr[7] = do_iso[shape];
	index_arr[8] = 0; // 0 = double, 1 = float
	index_arr[9] = 0; // Left empty for some future purpose
	for (lii1 = 0; lii1 < kintvars.num_li; lii1++)
	  index_arr[10+lii1] = tfvars.l_arr[kintvars.li_arr[lii1]];
	dimb = count * np*np*np * ni*ni*(ni+1)/2;
	bispec = (double*)malloc(dimb *sizeof(double));
	if (bispec == NULL){
	  fprintf(stderr, "Memory allocation failed.\n");
	  exit(1);
	}
	index = 10 + kintvars.num_li;
	count = 0;
	for (lii1 = 0; lii1 < kintvars.num_li; lii1++){
	  params.li[0] = lii1;
	  l[0] = tfvars.l_arr[kintvars.li_arr[lii1]];
	  for (lii2 = lii1; lii2 < kintvars.num_li; lii2++){
	    params.li[1] = lii2;
	    l[1] = tfvars.l_arr[kintvars.li_arr[lii2]];
	    index_arr[index] = count;
	    first = 1;
	    for (lii3 = lii2; lii3 < kintvars.num_li; lii3++)
	      if (triplet_relevant_wrap(TRIPLET_TYPE, kintvars.li_arr[lii1],
				    kintvars.li_arr[lii2],kintvars.li_arr[lii3],
				    stepstep, stepmax)){
		if (first){
		  li3min = lii3;
		  first = 0;
		}
		params.li[2] = lii3;
		l[2] = tfvars.l_arr[kintvars.li_arr[lii3]];
		switch (shape){
		case 0:
		  if (REALRUN==2)
		    bispectemp = compute_bispectrum_local(&params, shape);
		  if (REALRUN==3)
		    bispectemp = compute_bispectrum_local_indep(&params, shape);
		  break;
		case 1:
		  bispectemp = compute_bispectrum_equil(&params, shape);
		  break;
		case 2:
		  bispectemp = compute_bispectrum_ortho(&params, shape);
		  break;
		default:
		  fprintf(stderr, "Unknown shape: %d\n", shape);
		  exit(1);
		}
		for (i = 0; i < np*np*np*ni*ni*(ni+1)/2; i++){
		  bispec[i*index_arr[3]+count] = NORM*NORM * bispectemp[i];
		  bispec[i*index_arr[3]+count] *= ((double)(l[0]*(l[0]+1))
		    * (double)(l[1]*(l[1]+1)) * (double)(l[2]*(l[2]+1)))
		    / (l[0]*(l[0]+1) + l[1]*(l[1]+1) + l[2]*(l[2]+1));
		}
		free(bispectemp);
		count++;
	      }
	    index_arr[index+1] = kintvars.li_arr[li3min];
	    index_arr[index+2] = kintvars.li_arr[params.li[2]];
	    index += 3;
	  }
	}
	fp = fopen(filename_index, "wb");
	fwrite(index_arr, sizeof(int), dimi, fp);
	fclose(fp);
	fp = fopen(filename_out, "wb");
	fwrite(bispec, sizeof(double), dimb, fp);
	fclose(fp);
	printf(" *** Successfully create bispectrum. *** \n");
	printf("ni = %d\n", ni);
      }

  free_all_splines();
  free_tfvars();
  free(kintvars.li_arr);
  free(kintvars.x_arr);
  free(kintvars.w_arr);
  free(kintvars.alpha_arr);
  free(kintvars.beta_adi_arr);
  free(kintvars.beta_iso_arr);
  if (do_shape > 1){
    free(kintvars.gamma_arr);
    free(kintvars.delta_arr);
  }
  free(cl);
  free(bispec);
  free(index_arr);
  free(filenames_tf);
  free(iso_string);
  free(pol_string);
  return 0;
}
